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We determine the trajectories of vortex singularities that arise after a single vortex is broken by a 
discretely symmetric impulse in the context of Bose-Einstein condensates in a harmonic trap. The 
dynamics of these singularities are analyzed to determine the form of the imprinted motion. We 
find that the symmetry-breaking process introduces two effective forces: a repulsive harmonic force 
that causes the daughter trajectories to be ejected from the parent singularity, and a Magnus force 
that introduces a torque about the a^xis of symmetry. For the analytical non-interacting case we 
find that the parent singularity is reconstructed from the daughter singularities after one period 
of the trapping frequency. The interactions between singularities in the weakly interacting system 
do not allow the parent vortex to be reconstructed. Analytic trajectories were compared to the 
actual minima of the wavefunction, showing less than 0.5% error for small impulse strengths over 
the entire trajectory. We show that these solutions are valid within the impulse regime for various 
impulse strengths using numerical integration of the Gross-Pitaevskii equation. We also show that 
the actual duration of the symmetry breaking potential does not significantly change the dynamics 
of the system as long as the impulse strength is small. 



I. INTRODUCTION 

Vortices are ubiquitous to many diverse branches of 
science, such as fluid dynamics l[, meteorology [1], cos- 
mology d, Hquid crystals [5|, superconductivity 
II], solid state physics 0, and nonlinear singular op- 
tics [13, El- Vortices have been realized experimen- 
tally in Bose-Einstein condensates (BEG), obtained when 
bosons are cooled down to almost zero temperatures 
l22j . These vortices are expected to ofi^er interesting ap- 
plications in intcrferometry (23j , and as a means to study 
the behavior of random polynomial roots [l^]. Vortices 
are characterized by the presence of a phase singular- 
ity to which an integer number can be associated, called 
vorticity, topological charge, or winding number [25l.l26j. 
This singularity behaves as an individual physical entity 
whose motion can be studied separate from the conden- 
sate. The determination of this motion and its control 
is applicable in the study of many of the fields described 
above. Here, we determine the dynamics of these singu- 
larities when a highly charged vortex in a Bose-Einstein 
condensate is struck by a symmetry-breaking impulse. 

Specifically, we consider a dynamical situation in which 
a highly charged two-dimensional vortex is generated in 
an axisymmetric harmonic potential, and a symmetry- 
breaking potential is turned on transversely, i.e., in the 
plane containing the vortex, for a very short period of 
time, such that it can be described by an impulse with 
a potential V over a period of time At. This potential 
shows some rotational discrete point symmetry of order 
iV, that is, it reproduces itself under multiple integer ro- 
tations of 27r/iV • The topological charge of the vortex 
will experience a transformation, as discussed in (28l - l30| . 
We will show that the highly charged parent singularity 
will disintegrate into a number of single-charged daugh- 
ter singularities of different sign. The number and sign 
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FIG. 1: (Color online). Schematic of transformation pro- 
cess. The upper panels represent the potential along time 
while the lower ones are the corresponding phase singularity 
structure. The harmonic trap is represented in blue and the 
symmetry-breaking impulse with a green surface. The initial 
highly charged vortex, represented in e) as a red surface, is 
broken in five daughter singularities, one at the origin with 
charge —1 (blue surface) and four off-axis of charge -1-1 (red 
surfaces). The four off-axis singularities move outwards, even- 
tually coming back close to the origin due to the trap [black 
arrows in g) and h)]. 



of these daughter singularities are related to the pecu- 
liarities of the symmetry-breaking impulse |3l| . Figure [1] 
shows a representation of this discretely symmetric im- 
pulse acting on the parent singularity. 

Here we determine the path followed in a harmonic 
potential by the ejected daughter singularities after the 
impulse. We obtain these trajectories analytically for the 
non-interacting case by utilizing the Feynman propaga- 
tor for a harmonic potential, and determine their valid- 
ity in the weakly interacting case. In the non-interacting 
case we find that the parent singularity reconstructs it- 
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self from the daughter singularities after a period of time, 
i.e., there is a full quantum revival. In the interacting 
case quantum revival is blocked in the mean field the- 
ory: repulsion prevents the highly charged parent from 
being reconstructed, hence describing a helical trajectory 
around the origin. We discover that there are two effec- 
tive forces the singularities experience during symmetry 
breaking; a repulsive harmonic force that causes the tra- 
jectories to propagate outward, and a Magnus force that 
introduces a torque about the axis of symmetry. Our re- 
sults pave the way to the control and manipulation of the 
motion of singularities by means of symmetry-breaking 
impulses. The results are equally applicable to the neigh- 
boring field of nonlinear singular optics fllj by exchang- 
ing time evolution with axial-spatial evolution and the 
symmetry-breaking impulse with a inhomogeneous thin 
diffracting element. 

The study of dynamics of singularities and their inter- 
action is an exciting field with many potential applica- 
tions. The dynamics of vortex dipolcs; their interaction, 
oscillation, tunneling, and their collapse; has been theo- 
retically studied in the framework of BEC [s^ - lioj . Other 
structures of singularities and the interactions among 
them lead to elaborated trajectories [4l|, |42j , as discussed 
numerically in |4^. The geometry of vortex trajectories, 
like loops or hyperbolas, is related to vortex creation and 
annihilation and vortex interactions, and its study leads 
to a variety of vortex structures [25|, |43 - |47| . Moreover, 
the interpretation of the role of a phase singularity in 
quantum dynamics is an interesting issue, as well as the 
effect of th e dy namics of the singularities in the quantum 
system [1^,14^1 . Also, vortices in BECs, called vortex soli- 
tons in nonlinear optics, can show more than one off-axis 
singularity; in the latter case they are called vortex clus- 
ters. [26l. [50l - [56j . These structures are typically unstable, 
showing very slow dynamical decay rates, though some 
controversy has been built up around this issue [57l - l59| . 
Here, we obtain a breathing cluster of vortices, or plan- 
etary vortices, and we obtain numerical stability of this 
structure along the times of the evolution. 

In Sec. ini we define the system and the impulse, and 
briefly explain the transformation rule used to determine 
the symmetry of the post-impulse singularity structure. 
In Sec. Illli we introduce the initial highly charged vor- 
tex, with winding number £ = 3, that will be used in 
this research as our main test case, for brevity; similar 
behavior is seen for other initial winding numbers. This 
section also carries out the harmonic oscillator propa- 
gation integral used to determine the wavefunction after 
the impulse. Section HVl uses the wavefunction after sym- 
metry breaking to analytically derive the trajectories for 
the off-axis singularities that are broken out of the initial 
highly charged vortex by the impulse. Section fVl utilizes 
the trajectories to arrive at analytic descriptions of the 
equations of motion. These equations of motion are an- 
alyzed to understand the fundamental motion that the 
vortices undergo once symmetry is broken. Section IVII 
includes a comparison of the analytic trajectories with 



the local minima of the wavefunction. The analytic tra- 
jectories are also compared to numerical analysis of the 
Gross-Pitaevskii equation for various values of the non- 
linearity, g. In Sec. IVIIl we conclude. 

II. VORTICES IN A WEAKLY INTERACTING 
BOSE-EINSTEIN CONDENSATE 

Let us consider a system of weakly interacting bosons 
of mass M confined in a harmonic trap and condensed 
in the ground state at T <C Tbec, thus forming a BEC. 
We assume that one of the trapping frequencies is suffi- 
ciently high to reduce the dimensionality of the system 
to only two dimensions, but not near any potential res- 
onances. This system can be described by the Gross- 
Pitaevskii equation (GPE) 

ihdtij{Si, t) ^ Hi^iSi, t), (1) 

with 

tf = '^^' + Vi^^t)+93Bm'^ (2) 

where x = {x, y) G M?, and g^i^ is the coupling constant, 
or nonlinearity, defined by g^D = Airh^as/M , where Os 
is the scattering length and M is the reduced mass, and 
related to the effective interactions among the bosons in 
the trap. To model the symmetry breaking impulse, we 
consider a time-dependent potential given by 

r yo(x) Q<t<to 
t) = I yo(x) + AK)(x) <o < t < = + At 
{Vi{i) h<t. 

(3) 

We represent this potential in the upper panel of Fig. [T] 
We assume the length of the second region to be small. 
At <C 1, to model a Dirac delta impulse with constant 
area but short duration. Also, we assume that AVq is 
invariant under the action of the elements of a discrete 
rotational group Cn in order to view the effects of a 
symmetry-breaking impulse. The original and final me- 
dia own perfect rotational symmetry, denoted 0(2), with 
a potential given by 

Vo = ^ML,^{i^ +f). (4) 

Mathematically, we express the invariance property of 
the impulse as 

Ayo(Gx) = AK,(X) VGeCw. (5) 

Let us consider that for t < to the atoms arc condensed in 
a vortex of vorticity v, where v = 1 /2'k f-^ V0 • dl where T 
is a closed path encircling the axis of cylindrical symme- 
try of the vortex [2^. We assume that the phase singu- 
larity located at this axis is highly charged, v > 2. This 
vortex shows angular momentum £ equal to w, which is 
conserved if the symmetry is not externally broken [26j . 
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To analyze the effect of discrete symmetry poten- 
tials of order N in the properties of vortices, a quan- 
tity called angular pseudomomentum was associated to 
them [i^. To define this quantity, it should be noticed 
that any vortex solution of the GPE can be written as 
'tP{r,e) = e™'^u(r,6'), where u{r,9) = u{r,9+ ^), with 
m an integer. The effect of such a rotation is the adding 
of m times the same angle to its phase [s^, ■ It was 
shown that the values of m are restricted by the order of 
symmetry [2^ . and hence: 



0,±1,±2,...,4 



0,±1,±2,. 



jV-l 
2 



even N 



odd TV. 



(6) 



Also, it was shown that every vortex presents a phase 
singularity of charge m in the origin, where the charge of 
a singularity is Vj = 1/2ti V0 ■ d\ where Vj is a closed 

path that encircles only this singularity [2^. It was also 
shown that one can relate the angular momentum £ of 
a circularly symmetric vortex struck by a discretly sym- 
metric potential of order with the angular pseudomo- 
mentum m of the wave function in the discrete symmetry 
media by the transformation rule [28l - [30| : 



i-m = kN, 



(7) 



where k is an integer. The process of transformation has 
been described microscopically as the disintegration of 
the highly charged vortex into a number of smaller vor- 



tices |3l| . In the axis of symmetry, a vortex of charge m 
remains, according to [26j . while the integer k in Eq. ([7]) 
is related to the number of rings of single charged vor- 
tices emerging from the axis j3ll |. For example, consider 
a circularly symmetric vortex with vorticity w = 3 being 
broken by an iV = 4 symmetric impulse, the transfor- 
mation rule in Eq. ([7]), gives m = — 1 with fc = 1, and 
the central singularity has charge —1. In this case, the 
initial vortex will result in 5 singularities after symme- 
try breaking. One stays at the origin, with a new charge 
of = —1, and one ring of four Vj — +1 charged sin- 
gularities comes symmetrically off-axis, as represented in 
Fig. [TJ Here we are interested in the trajectories fol- 
lowed by these ejected daughter vortices after emerging 
from the axis. For illustration purposes, in the following 
we will consider only this particular case, even though 
all results can be extended to an arbitrary charge of the 
initial vortex and symmetry order of the impulse. 



III. PROPAGATION OF A VORTEX AFTER A 
SYMMETRY BREAKING IMPULSE 

In this section as well as Sec. II VI we consider the linear 
case, g = 0, and we use the substitutions 



(jjt — T, L = \ , X = Si/L, 

nvjj 



(8) 



and we take = 0. We consider the initial normalized 
vortex of charge ^ = 3 given by 



4>n, m(x) 



1/ . ^3 1 

6 V2"+™n!TO!7r 



X exp 



V 



H„{x)Hm.{y), 



(9) 



which is a solution of Eq. ([T]). In the following we consider 
the vortex with charge £ = 3 with less energy, i.e., we set 
n = m = 0. As shown in Appendix 1X1 the amplitude of 
the vortex wavefunction after the action of the impulse 
will be given by 



(10) 



where (/) ( represents the vortex wavefunction before the 
impulse [6l| and we use the complex notation x = x + iy, 
and At is the duration of the impulse. Using the Taylor 
expansion of the evolution operator given in Eq. (|B2p of 
Appendix IB] gives us 



_ iAT(i;ox*+uix"') 



0(X). 



(11) 



If we carry out another Taylor expansion for the expo- 
nential, we get the expression for (f>(x) that will be used 
for the propagation in the final medium 

0(X) = (1 + ivoArx^ + iviATx*'')<f>ix)- (12) 

Inserting the initial state 0„„j given in Eq. Q into the 
previous equation, we arrive at the final form of our wave- 
function after symmetry breaking. 

0(x) = [l + ivQ{x + iyY + ivi{x - iyY] 
^/ • ^3 1 

6 V2"+™n!TO!7r 



X exp 



y 



Hn{x)H„,{y), (13) 



where we have absorbed the length of the impulse, Ar, 
into the constants vq and vi such that vqAt = vq and 
171 At = vi- The two-dimensional harmonic oscillator 
propagator given by 



^(x, t) 



1 



V'(a;o, T-o)exp 



2i7rsinr 

i 



exp 



t cos Tlx 



dxodyo 



2 sin T 



2 siuT 

[i^l + yl) cos T - 2{xxq + yyo)] 

(14) 



valid for t > 0, is then used to evolve the initial func- 
tion Eq. ()13p in the transverse plane. Taking as the ini- 
tial state 000 J we get for the vortex field after symmetry 
breaking 



0c(x, t) = e 



^Ut-\xVI2 



(15) 



+ ivrx* [-24 + 246^'- + \x\\\xY - 6)' 
+36e4-(|xP - 2) + 12e2- {6 + |xP(|xP " 6)}] } 
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FIG. 2: The wave function (pix, y, r) for different times, (a) 
to ((f) amplitude and a corresponding zoom for different time 
slices, (g) to (i) corresponding amplitude. Red dot represents 
position of negatively charged singularity at the origin, yellow 
dots represent the positions of the single positively charged 
singularities. 



where we have used the identities x^/x* = X^/lxP ^^^d 
X^/x* = x''/lxPj E^nd where 



Fix, r) = 



6^ 



(19) 



It becomes immediately apparent that F(x, r) is C4 
invariant due to the dependence on only x'^ and terms. 
Because F{x, t) is C4 invariant, 



'/'(ex, t) 



V(X, t). 



(20) 



where e 



is the elementary rotation of 4th order. 



Thus, as expected from the analysis in Ref. [26|, [291, |30|, 
l60j . using the transformation rule, the solution preserves 
the winding number m — —1 for the center singularity. 

In Fig. [2] (a), (b), and (c) we represent the amplitude 
of this function for different times. A closer view of these 
amplitudes is shown in Fig. [5] (d), (c), and (f), while 
Fig-[l](g)> (h), and (i) show the phase for the same times 
and also near the origin. In this last figures the position 
of the singularities is highlighted. The positively charged 
singularities are ejected from the origin and follow some 
trajectory. In Sec. IIV[ we find the expressions for these 
trajectories. In Sec. |V] we will determine the effective 
forces acting on these singularities. 



IV. TRAJECTORIES OF THE PHASE 
SINGULARITIES 



Let us write Eq. (|T5|) as 

0c(x, r) =e-«— 



6 



(16) 



(Ao(t)x^+A+x' + A_(|x|, r)x*), 



where 



and 



4ir 



A_(|x|, r)=ii,i(-24 + 24e«- + |xP(|xP-6)2 (17) 
+ 366^^^(1x1' - 2) + 12e2-(6 + |x|'(|xl' - 6))). 



The expression in ([T6|) has the form predicted by our 
previous symmetry arguments (26j since it can be written 

as 



g-8ir-|x|72^* 



A+x' + Ao{t)x' 



+ ^-(1x1, r) 



X*F{x, r). 



(18) 



Let us obtain the trajectories followed by the phase 
singularities by finding the zeros of the complex wave- 
function (j){x, t). From ([T5)) . we see that there are two 
situations when 4'{x, t) ~ 0, when X* = and when 
Fix, t)^0. ^ 

For the former situation, we study the behavior of the 
wavefunction near the origin by developing 4>{x, ''") in a 
Taylor series around x = 0, obtaining: 



0(X, r) « 32e- 



-5iT 



(21) 



Evidently, the singularity at the origin is due to the sym- 
metry breaking of the initial vortex into C4, as seen by 
the dependence on the symmetry breaking parameter vi . 
We see again that this singularity has winding number 
TO = —1, as evidenced by the factor of x*- If we set 
fi = 0, then ^_(|x|, t) = 0, and the expansion about 
X = is instead 



Hx, 



■'^^x^ 



(22) 



which preserves the initial winding number of Z = 3, as 
seen by x^- 

The latter type of phase singularity, when F{x, r) = 0, 
is more difficult to analyze because we have to work with 
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the complex roots of the nonhnear equation -F(x, ■'') = 0. 
This is the same as solving the equation 



(23) 



To make the calculation easier, we assume that the two 
symmetry breaking parameters are equal. Thus, we take 

Vo=vi= V. 

If we go to the v = limit, we see that = and 
t) = 0. For F{x: r) = to be true in this limit, 
it follows that as v — 0, Ao(''")lxl'* 0, and therefore 
X — >■ 0, leading to the conclusion that x = x{'^)j a-nd the 
statement that in the u <C 1 regime, |x| ^ 1- 

For small values of Xi the first terms that reappear in 
([23]) are those in A_(|x|, t) that depend on |xP- Due to 
Ixl being much less than 1, it follows that |xP > Ixl'' > 
Ixl^- By expanding out yl_(|x|, t), we see that the |xP 
term is 



lim A_(|x|, t) w \xW-24iv + 72ie^"v 



(24) 



Using this approximation, we can instead solve the equa- 
tion 



^o(t)x' + ^-„«i(|x|, t) =0, 



(25) 



where we have kept only the nonzero terms from the v <C 
1 limit. 

Thus, to order |xP, 



e^'^X^ + \xW-24iv + 72ie^'^v 
~ 72ie'^"v + 2'^ie^"v) = 0. 

If we solve for x^, 



(26) 



24w - 72ie2*^z; + 72ie'**^z; - 24^6^*^^ ^ 



up(t)|x|^ 



(27) 



The simplest mathematical object to calculate now is |x|- 
This is done by taking the modulus of the previous ex- 
pression and dividing by |xp. We obtain 



IXl 



192i)i sin^(r). 



(28) 



The equation above provides the evolution of the ra- 
dial coordinate of the off-axis phase singularities after 
the symmetry-breaking impulse has been applied. Recall 
that |x(''')P = x{t)'^ + viT)^ = ^{t)'^ so that in polar 
coordinates the radius of the phase singularity trajectory 
is given by 



R{T)^8Vsivi sin3(T))i/2. 



(29) 



To find 0{t), we need to look back at (pTj) . If we rewrite 
X and p in modulus-argument complex form, x becomes 
|x|e'^ and p becomes |p|e*'''. Equation (P7)) becomes 



\x\'e 



2 i4e 



v\p\e 



(30) 
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FIG. 3: Trajectories followed by ejected daughter singularities 
after symmetry breaking. 



We saw in (P?]) that |xP = vp, so the previous equation 
becomes 



j4e 



= 46' = 7 + 2mT. 



(31) 



Thus, the evolution of the polar coordinates of the phase 
singularities is provided by the phase of p(t). From (|27p . 



p(r) 



24i - 72ie2'^ + 72ie*''' - 24ie^ 



(32) 



The phase of p{t) is found by taking the arctangent of 
p(t). This gives us 



2n7r -|- arctan 



sin*(T) 



— cos(t) sin (r) 



, , mr T 

(r) « . 

^ ' 2 4 



(33) 



The analytic trajectories calculated above are repre- 
sented in Fig. [21 



V. EQUATIONS OF MOTION 

In this section, we analyze the trajectories found pre- 
viously to determine the effective forces acting on the 
singularities during and after symmetry-breaking. The 
equation of motion for the complex coordinate i?(r)e'^'^'^^ 
shows that the system undergoes nontrivial motion corre- 
sponding to a harmonic oscillator with complex and time- 
dependent frequency. This leads us the discovery of an 
effective singular force that causes the daughter singular- 
ities to initially repel. Due to the complex nature of this 
effective potential, we show that the system of vortices 
is non-conservative. We then describe the torque about 
the origin acquired by the singularities during symmetry- 
breaking. This torque arises as the product of the Mag- 
nus force. 
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A. Radial and Angular Equations of Motion 

Now that we have expressions for R{t) and 0{t), we 
can find the equations of motion for the off-axis singular- 
ities. By taking the derivatives of R{t), we see that the 
velocity and acceleration in the radial component can be 
expressed as 

i?(r) =12\/3cos(r)(i;isin(r))i/2 (34) 
R{t) =3%/3(3 cos(2t) - csc(r))i/^ (35) 

By taking the derivatives of 0{t), we see that the angu- 
lar velocity is constant, 9{t) = —j, and therefore the 
angular acceleration is zero. 

However, the fact that there is an angular velocity at 
all tells us that the singularities have acquired a torque 
about the axis at some point in their creation and prop- 
agation. If we recombine R{t) and 0{t) into the complex 
coordinate xi''') = R(T)e^^''^'> once again, we can study 
the behavior of the singularities immediately after sym- 
metry breaking. 

If we Taylor expand equations ([^ and ([55)1 around 
r = 0, i.e immediately after symmetry breaking, we get 



R{t) « 8V3v^T3/^ 



and 



(t) 



TT 
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(36) 



(37) 



The previous expansions give us the complex coordinate 
x(r) right after symmetry breaking, such that 



X(r) 



(38) 



We now proceed to derive the equation of motion asso- 
ciated to (|55)) . If we differentiate (pS)) with respect to r, 
we see that 



Xir) = ( I - i ) Xir) 



(39) 



A second derivative of psp will provide us with the equa- 
tion of motion in complex notation: 



12-t(12z + t) 
X (j) = T7n? xir). 



16r2 



(40) 



If we let = i 



as 



and r2? = -j^, we can rewrite 



x"{r) = -{nl + ^n\)x{T). 



(41) 



Evidently, the phase singularities experience a nontrivial 
type of force. The previous equation represents a spe- 
cial type of harmonic oscillator in which the frequency is 
both complex and time dependent. Since the frequency 
is complex, we do not expect the system to be conserva- 
tive. We can prove this statement by manipulating PO)) 
and its conjugate in the same manner we would do to 



establish conservation of energy in a standard harmonic 
oscillator. First, we write the conjugate of ([Ml : 



(42) 



Next, we muhiply (gO]) by x'*{'r) and ^ by x'(t) and 
add the two resulting equations to obtain 

X"(r)x'*(r)+X'(r)x"*(r) = 

= -{nl + inl)x{T)x'*{T) + ^{T){-nl + i^I\)x*{t) 
= -17§x(T)x'*(r) - ^n\x{T)x'*{T) nlx'ir)x*ir) 

+ inlx'{r)x*{r) 
= -nUx{r)x'*{r)+x'{r)x*{r))+^nl{x'{r)x*{r) 

-x(t)x'*(t)). 

We immediately recognize that the left hand side and 
the first term of the right hand side are total derivatives. 
If we rewrite the total derivatives, we get 



^{X'{r)x'*{r)) = -nl^ix{r)x*ir)) 



dr 



d_ 



+ zfif(x'(r)x*(r)-x(r)x'*M). (43) 



Unfortunately, fig is time dependent, so we cannot just 
combine the total derivatives. Instead, we must subtract 
the term with -§p^o- The total derivative of the fip term 
is 

^inorxiT)x*ir)) = x{r)x*{r)^nl+nl^x{r)x*ir). 

(44) 

This allows us to rewrite as 

^ (X'(r)x'* (r)) + ixir)x* (r)) (45) 

dr dr 

= zl7?(x'(r)x*(r)-x(r)x'*(r)), 



^ (X'(r)x'*(r) + nlxir)x*{r)) x(r)x*(r)^f^g 
= »f^?(x'(r)x*(r)-x(r)x'*(r)). 

If we replace x'i^) and x'*(j') with their functional val- 
ues, and evaluate the derivative of f2§, we obtain 



^{x{r)x'*iT)+nlxir)x*ir))^ 

V2r3 2 



(46) 
X(r)x*(r). 



If we define the energy of the system the same way we 
would a typical harmonic oscillator, 

E = ^x'ir)x'*ir) + lnlxir)x*ir), (47) 

it is clear that there is gain in the system. We can see 
the value of the gain by considering the derivative of the 
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energy: 



dE_ 
d7 



2dT 



-J^ + ^]xir)x*{r) 



At 



1 



1 



lx(r)|'>0. 



(48) 
(49) 

(50) 



Thus, energy is not conserved by our equations of mo- 
tion governing singularity or vortex motion. However, 
the GPE does conserve energy. Thus energy is being ex- 
changed between the singularities and the remainder of 
the Bose-Einstein condensate described by the full GPE. 

The presence of this effective harmonic motion explains 
why the post-symmetry breaking singularities expel from 
the origin. However, we still need to explain the effective 
torque that the singularities seem to experience. To un- 
derstand this torque better, let us rewrite our complex 
coordinate in Cartesian coordinates via the definition of 
x(t) = x{t) + iyij). This gives us 



x"{r)^x"{T)+iy"{T) 

= -{nl + inl){x{T\y{T)) 

= -{^il + in{)x{T) 

= -nlx{T) - m\x{T) - n - 



(51) 



If we collect the real and imaginary parts, we arrive at 

x"{t) = -nlx{T) + nly{T) (52) 
y"iT) = -nly{T)-nixiT). (53) 

We can write the previous equations in vector form as 



r"(T) = -02r(r)+02 



1 

-1 



r(r). 



(54) 



In order to see how a torque comes into our system, we 
need to rewrite the ftf matrix term in three-dimensions 
(3D). To do this, we construct the external 3D vector 
A = (0, 0, nj) such that 



r X A 



1 J 

X y 



ni 



nl{y, -X, 0) = nf 





-1 



(55) 

where VTir) is the transverse plane and i, j, k are 
unit vectors, which we are working in. Therefore, the 
equation of motion for phase singularities can be rep- 
resented in 3D, although the motion is restricted to a 
two-dimensional plane r(r) — (x, y, 0)[65j. We write our 
3D representation as 



(56) 



This equation of motion shows the simultaneous presence 
of a harmonic force and an external force associated with 
a torque. The fact that the latter is associated with a 



torque can be checked by calculating its effect on the 
angular momentum of the phase singularity - L = r x r'. 
If we look at the derivative of the angular momentum, 
we see that 



dL _ d 
dr dr 



(r X r') 



r X r + r X r 



(57) 



If we evaluate this cross product using our expression for 
r"(T) in dSSl), 



r X r" = r X (-ilgr + (r x A) 

= -^l{r X r) + r X (r X A) 

= rx(rxA). (58) 

Using the vector triple product [62j . we obtain 

r X r" = r(r • A) - A(r • r). (59) 

Because A is defined only to have a z component, and 
our position vector is two-dimensional, the dot product 
of r with A vanishes, leaving 

r X r" = — A(r • r) 
= -A|rp 



(0, o,~nl\rTf). 



(60) 



Finally, we arrive at 



dL 

dT 



r= (0, 0, -nllrrl 



0, 0, 



4t 



(61) 

The previous equation shows that the angular momen- 
tum has variance only in the z direction, which means 
there is a torque that causes rotation in the x, y plane, 
as we expect. Because the value of the torque is negative, 
our singularities rotate about the origin in a clockwise 
manner, as our trajectories in Sec. IIVI were seen to do in 
Fig. O In Appendix [Cl we check that the calculation of 
the energy is correct using our 3D formalism. 

Let us note that all of these results apply, when prop- 
erly rotated, to any of the four phase singularities mov- 
ing away from the center of symmetry. This is due to the 
four-fold symmetry of our solutions and it is reflected in 
the four solutions that we have for the angular coordinate 

e{T) in (1331). 



B. 



Dynamics of Phase Singularities in Free Space 
after Rotational Symmetry Breaking 



In Sec. IV Al we derived the equation of motion for the 
four phase singularities that arise immediately after sym- 
metry breaking by a discretely symmetric impulse. We 
found that the breaking of rotational symmetry causes a 
vortex to cluster in a central singularity carrying topo- 
logical charge equal to the angular pseudo-momentum m 
and a "wave" of N (N being the order of symmetry of the 
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impulse) single charged phase singularities with particle- 
like motion moving away from the symmetry axis. The 
dynamics of these phase singularities as point-like parti- 
cles is described by the equation of motion in ([55)) (for 
the case iV = 4). This equation is very interesting be- 
cause it shows that, despite the wave function describ- 
ing the propagation of matter corresponds to linear har- 
monic propagation, the clustered phase singularities do 
not move as harmonic oscillators. In fact, right after the 
action of the impulse, they experience two types of forces, 
as described by the right-hand side of ([55]): 



A harmonic repulsive force given by il^r^T). 

A rotational force F ~ (r(T) x A) generating a 



torque M 



^A. 



Both forces have a peculiar behavior, 
them separately. 



Let us analyze 



1. Effective Harmonic Potential 

This effective potential is crucial because it is responsi- 
ble for the dissociation of the initial highly-charged vor- 
tex with topological charge £ = 3 into the central sin- 
gularity of charge 7ti = — 1 and four vortices of charge 
Vj = If the interaction was attractive, the four vor- 
tices would remain at the origin (the center of symmetry) 
since both the initial position and initial velocity are zero. 
However, we find that the interaction is repulsive because 
rig = — ^ < for small values of t. Nevertheless, 
a repulsive harmonic interaction is not enough to guar- 
antee the motion of the broken singularities away from 
the origin since their position and velocity are initially 
zero. They would remain there in a situation of unsta- 
ble equilibrium since the force upon them would be zero. 
Something else is needed to trigger the expansive mo- 
tion of the broken singularities. The mechanism is the 
existence of a nonzero, in this case singular, repulsive 
potential at t = 0. 



'H 







(l6 


4t2 J 



r(r) ^ -^T 



3/2 



1 



oo. 



(62) 

If we analyze the form of the effective harmonic potential 
for small values of i, we see from (|48l) that 



VH{r) 



1 



3 



rr i < 1, 



(63) 



indicating the presence of a singular repulsive potential 
at t = 0. The curvature of the quadratic potential is, 
thus, infinite and negative right after the symmetry is 
broken, so the force on the escaping singularities is non- 
zero when they are located at the origin when t = 0. This 
singular potential is the reason why the singularities start 
to move away from the center of symmetry. The fact 
that the potential and force are singular at t = does 
not produce any issues in the velocity and position of the 



fleeing singularities when t = because the acceleration, 
which has the form r"(r) ^ has first and second 

integrals of the form: 

r3/2 + C', 



r'(T) ^ \/t + and r(r) 



(64) 



which arc both finite at t = and compatible with the 
initial condition r'(0) = and 7'(0) = when the con- 
stants are taken to be zero. 



2. Torque 

As seen in (|61[) , there is an r dependence in the torque 
that the singularities experience around the origin once 
symmetry is broken. The torque is zero when t ~ since 
T ~ izt2^ Q due to the initial condition of r(0) = 0. 
Thus, the singularities must start moving away from each 
other, making r ^ 0, before the external torque can take 
effect. This allows us to conclude that the singular repul- 
sive effective harmonic potential acts on the singularities 
before they can acquire any angular momentum. 

As our vortices acquire a linear velocity away from the 
origin, they become subject to the Magnus effect. This 
effect creates a force perpendicular to the direction of 
motion according to -F = Sicj x i/), where 5 is a property 
of the medium the vortex is traveling through. oJ is the 
angular rotational velocity of the spinning object, and v 
is the linear velocity. This perpendicular Magnus force 
causes the vortices to follow a curved path. If we evaluate 
the expression for the Magnus force, we see that 

Siijx v) 
= S{0, 0, Lj) X (^x, 



(65) 



where the velocity vector was formed by taking the 
derivative of the position at small t given in (|36p . If we 
evaluate the torque associated with the Magnus force, we 
see that 

T — fx F 



{x, y, 0)x S [ -^y, ^2;,0 



0, 0, 



35a; 2 



(66) 



which is consistent with our expression for the torque 
found in the previous section with S = —I and u = 1, 
verifying that the torque associated with the singulari- 
ties after symmetry breaking is generated by the Magnus 
force. 

Therefore, the dynamics of our singularities after sym- 
metry breaking can be described as follows: First, the ac- 
tion of the symmetric impulse introduces an effective sin- 
gular repulsive harmonic potential that splits iV = 4 sin- 
gle phase singularities out of the original highly charged 
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vortex. As these singularities begin to travel away from 
the origin, they gain angular momentum from the effec- 
tive external torque and rotate around the axis of sym- 
metry. Eventually, the effective harmonic potential is 
overpowered by the trapping potential, so the singulari- 
ties travel back toward the origin, and settle into oscilla- 
tory motion about the origin. Eventually, they fuse back 
to the origin and reconstruct the initial vortex for the 
non-interacting case only. 

We can also find an expression for the maximum radius 
the singularities achieve as a function of given by 

-Rmax - 8V3^, (67) 

confirming that the strength of the impulse directly af- 
fects the motion of the singularities. 

VI. TRAJECTORIES OF SINGULARITIES IN 
THE PRESENCE OF INTERACTIONS 



y 




X 



FIG. 4: (Color online). The calculated trajectories (pink) 
are plotted against the amplitude minima (blue) for discrete 
time steps, using a value of v = 0.00005 for the area of 
the symmetry-breaking impulse. The approximated trajecto- 
ries are a good description of the singularity motion for small 
impulse areas. 



We can compare the calculated trajectories from sec- 
tion |IV] with the actual minima of the wavefunction am- 
plitude as well as numerical data generated by numeri- 
cally solving the GPE ([2]) . To begin, we compare the an- 
alytic solutions with the actual minima of the wavefunc- 
tion to determine the accuracy of the calculated trajecto- 
ries for various impulse strengths. Later, we compare the 
analytic trajectories, solved for a non-interacting BEC, 
and compare them to the numerical data for the same 
symmetry-breaking process in a weakly interacting BEC 
with various particle interaction strengths to determine 
the validity of our results in the nonlinear case. 



A. Linear Comparison 

To compare the analytic trajectories to their actual lo- 
cations in the wavefunction, we must find a way to track 
the singularities. Due to the non-analyticity of the wave- 
function, we must use the Minimize command in Mathe- 
matica in order to track the singularities for various time 
steps. Using a Do loop, we can append the location of 
the minima in the fourth quadrant to a list and plot the 
trajectories. As we do so, we can compare the calculated 
trajectories (in pink) with the located minima (blue) for 
V = 0.00005 in Fig. |4l The value of w is a numerical 
representation of the impulse area, AVoAt, as described 
by the potential profile in [ITl 

We can calculate the error between the two trajectories 
by using the formula 



Log. 



10 



^calc ^theory 



2 (^calc 



rth, 



(68) 



where 



y 



As one can see in Fig. [Sj the error is 
largest at the apex of the petal loop, but is still within 
0.5% of the amplitude minima for v = 0.00005. The 
first few points in Fig. [S] have very large error due to 



the close proximity of the singularities immediately after 
symmetry breaking. The minimize command searches 
for a local minima, so when all four external singularities 
are very close to each other, the local minima could be 
from any of the singularities, increasing the error. 



0.005 




1 X 10-* 
5 X 10"' 



1 X 10"' 



FIG. 5: (Color online). We plot the error on a logarithmic 
scale between the calculated trajectories and the local minima 
of the wavefunction amplitude for v = 0.00005. The error is 
always less than 0.5%. The scattered points at the beginning 
and end of the plot are due to the minimize command in Math- 
ematica being unable to distinguish between the four external 
singularities when they are extremely close to the origin. 



The error is insignificant until the outer edge of the 
petals. The increase in error is most likely due to the 
various approximation techniques used to calculate the 
analytical trajectories, one of which was working close to 
the origin. 

We can increase the value of the impulse area, v, to 
observe the loss of validity as the duration of the impulse 
increases. As one can see in Fig. [BJ the error signifi- 
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cantly increases as the impulse area, v increases. The 
error becomes greatest near the apex of the petal struc- 
ture, while still being within 5% near the origin. The 
error rises above 5% for times greater than r « 0.5 for 
the large impulse area, v = 0.005. 





V(x,y)=0.5 
V(x,y)=0.05 
V(x,y)=0.005 
Analytic 



FIG. 7: (Color online). We plot the numerical trajectories for 
impulse durations of At = 0.1, At = 0.01, and At = 0.001, 
and compare to analytic trajectories. We see no significant 
change as the duration is increased. 



FIG. 6: (Color online). We plot the error for impulse areas 
of V = 0.005 (yellow), v = 0.0005 (pmk), and v = 0.00005 
(blue). The error increases significantly as the duration of ... . 

the impulse increases, but remains the lowest near the origin. 1.000- ' '•' 

0.500 



These results show that the analytic trajectories are 
the best approximations for very small impulse areas, but 
are still valid near the origin for larger impulse areas. 



B. Impulse Approximation 



O.IOO- 
0.050 



0.010 
0.005 



To determine whether it is the duration of the impulse 
or the height of the impulse that affects the error, we 
include numerical integration studies similar to those in 
Sec: I VI CI but for g = 0. Three studies were evaluated for 
impulse area v = 0.0005, each with a different duration 
and height. 

V{x, y) = 0.005 At = 0.1 v = 0.0005 (69) 
V{x, y) ^ 0.05 At = 0.01 v = 0.0005 (70) 
V{x, y) = 0.5 At = 0.001 v ^ 0.0005 (71) 

If we plot each case, we see that there is no significant 
difference as the duration of the impulse is increased, as 
shown in Fig. [T] 

If we calculate the error between the analytic and nu- 
merical trajectories, we see that the error stays below 
10% once the singularities leave the origin. The large er- 
ror near the origin is due to the closeness of all off-axis 
singularities to the central one, together with the im- 
possibility of locating the singularities with an accuracy 
smaller than the grid spacing used in the numerical simu- 
lations of the Eq. ^ . We plot the error on a logarithmic 
scale in Fig. [51 

These results allow us to conclude that the actual du- 
ration of the impulse does not significantly change the 



^ 20 40 60 80 Too 120 

FIG. 8: (Color online). We plot the error between numerical 
trajectories for impulse durations of At = 0.1, At = 0.01, and 
At = 0.001, and analytic trajectories. We see no significant 
change as the duration is increased. 



dynamics of the system so long as the total impulse area 
is small. For the analytic trajectories, this means the ap- 
proximation is valid to within 5% for impulse areas less 
than vAt = 0.0005. For small areas, the impulse only 
serves to break the symmetry of the singularities, and be- 
comes negligible if the duration is increased. In essence, a 
shallow potential for a longer time which is not governed 
by the impulse approximation has the same effect as an 
extremely strong potential for an infinitesimal amount of 
time, as long as the area vAt remains small. 



C. Numerical Comparison - Nonlinear Case 

The GPE, given by was solved numerically for the 
same impulse used in the analytical analysis for various 
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values of the nonlinearity, g. The nonlinearity depends 
explicitly on the scattering length between particles. The 
time at which the nonlinearity becomes significant is ap- 
proximately 



_ hL 
Ng 



(72) 



derived by units considerations from the renormal- 
izcd effective 2D interaction strengthjol] g = g2D = 
\/%TiWUlMas oc gw, with Uz the transverse harmonic 
oscillator frequency, and = yjh/raujz', we use simply 
"g" for our 2D effective interaction strength for simplicity 
of notation. We observe the numerical data for attractive 
nonlinearity and repulsive nonlinearity to see the struc- 
ture of the vortex trajectories in each case. The effect we 
pursue is in the very core of the vortex, but there is an 
unavoidable limitation related to the grid spacing neces- 
sary to compute the minima of the wavefunction. Thus, 
the numerics have large error near the origin where the 
vortex cores initially overlap and again approach closely 
at later times. 



1. Repulsive Nonlinearity 

Repulsive nonlinearity arises when the particles in a 
Bose-Einstein condensate interact with one another via 
a positive s-wave scattering length, corresponding to pos- 
itive values of g. As the nonlinearity becomes larger, the 
trajectories begin to interact at further distances from 
each other. This prevents the singularities from recom- 
bining at the origin, and instead sends the trajectories 
into repulsive motion, similar to the behavior of like- 
charged particles, before returning to the oscillatory path 
about the origin. Trajectories derived from numerical in- 
tegration for various repulsive nonlinearities, g, can be 
seen in Fig. El 

By increasing the nonlinearity from g = io g = 
1, 2, 3, 4, wc sec that when the singularities come back to 
the origin, the nonlinearity begins to show its effects, as 
seen by the paths taken by the numerical data. As seen in 
the previous figure, once nonlinearity is introduced, the 
singularities interact before traveling straight across the 
origin. In the non-interacting case, the singularities do 
not come back to the origin, but instead switch directly 
to another of the four loops. To show this, we plot in Fig. 
|9]the trajectories for all four singularities in linear case, 
while only the trajectory of one of them in the nonlinear 
ones. 



2. Attractive Nonlinearity 

Attractive nonlinearity corresponds to a negative value 
for the s-wave scattering length, resulting in negative val- 
ues of the nonlinearity, g. As the attractive nonlinear- 
ity becomes increasingly negative, the trajectories begin 
to interact at further distances from each other. This 
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FIG. 9: (Color online). Comparison between repulsive non- 
linear numerical data for g = 1, 2, 3, 4. As the nonlinearity 
increases, the singularities are less likely to return to the ori- 
gin. The repulsive nonlinearity sets the trajectories in a com- 
pletely different orbit than the nonlinear case of g = 0. This 
repulsive motion is similar to a system of like-charged parti- 
cles in that the singularities interact with each other before 
returning to the oscillatory path about the origin. 



prevents the singularities from recombining at the ori- 
gin, behaving similarly to a system of planets, essentially 
"sling-shotting" around each other before returning to 
their oscillatory paths. Trajectories calculated numeri- 
cally for attractive nonlinearity can be seen in Fig. [TUl 

By increasing the attractive nonlinearity from g = to 
g = —1, —2, —3, —4, we see that the nonlinearity begins 
to show its effects near the origin, as seen by the paths 
taken by the numerical data. The negative nonlinear- 
ity sets the trajectories into orbital motion, essentially 
"sling-shotting" around each other before returning to 
oscillatory motion. However, with the attractive non- 
linearity, the trajectories return to the opposite orbit as 
with the linear g = data. Again, in Fig. [10] we show 
the trajectories for all four singularities in the linear case, 
and only one of them in the nonlinear ones. 



VII. CONCLUSIONS 

We have analytically described the equations of motion 
for the off-axis singularities that arise after the action of 
a symmetry-breaking impulse on an initial single highly 
charged vortex. For an initial vortex of vorticity i = S 
at the origin and a C4 discretely symmetric impulse, the 
symmetry of the initial vortex is broken into C4 as well. 
Four vortices with charge Vj = oscillate about the 
origin in a flowering pattern. A single vortex of charge 
Vj = —1 remains stationary at the origin. 

All future evolution of the singularities is determined 
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FIG. 10: (Color online). Comparison between attractive non- 
linear numerical data for g = —1, —2, —3, —4. As the non- 
linearity increases, the singularities are less likely to return 
to the origin. The negative nonlinearity sets the trajectories 
into orbital motion in the opposite direction as the nonlinear 
case, essentially "sling-shotting" around each other, causing 
the singularities to interact with each other before returning 
to the oscillatory path about the origin. 



by the order of symmetry of the impulse. The singular- 
ities are imprinted by the impulse and "remember" the 
effect of symmetry-breaking once back to an ordinary 
confining harmonic potential. It is interesting to note 
that the actual form of the impulse does not change the 
motion of the singularities. It is the order of symmetry 
that determines all future propagation patterns. 

The calculated trajectories of the off-axis singularities 
give rise to a blossoming structure. The singularities pe- 



riodically oscillate about the origin, while rotating about 
the axis of symmetry. The disassociation of the initial 
highly-charged vortex into several smaller vortices is due 
to an effective singular repulsive harmonic potential that 
is introduced by the symmetry-breaking impulse. The 
singularities also acquire angular momentum around the 
axis of symmetry due to an external effective torque 
caused by the Magnus force. Once the effective repul- 
sive potential is overpowered by the trapping harmonic 
potential, the singularities settle into an oscillatory pat- 
tern as expected in a harmonic trap. 

The analytic trajectories were compared with the lo- 
cal minima of the wavefunction for impulse strengths 
of V = 0.005, V = 0.0005, and v = 0.00005. Com- 
parison with the local minima showed the trajectories 
to be within 0.5% error for v = 0.00005, 5% error for 
V = 0.0005, and 50% error for v = 0.005. The increase 
in error as the impulse duration is increased is due pri- 
marily to approximations made in the analytical analysis. 
By superimposing the analytic trajectories with the lo- 
cal minima, we see that they are in agreement for small 
impulse strengths. The actual duration of the impulse 
docs not significantly change the dynamics of the sys- 
tem so long as the total impulse area is small, less than 
V{x, y)At = 0.0005. 

The initial break-up of the singularity is completely 
controlled by linear effects. It is only long-time behavior 
that requires full nonlinear analysis due to the interaction 
between particles in an interacting BEC. 
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Appendix A: Evolution operator 

For the potential profile given in ^ , we can decompose 
the evolution operator into three separate operators, one 
for each region, according to 



-Ti)^iHAT^lHTo 



(Al) 



Now. H = Hho+Vsb, where Hho is the initial Hamilto- 
nian of the harmonic oscillator, and Vsb is the symmetry- 
breaking potential, both in dimensionless units. We de- 
fine Hq = Hho + Vo and Hi = Hho + Vi- According 
to the potential profile given in Eq. ([3]), the evolution 
operator can be rewritten as 



Let us analyze the evolution operator for the impulse. 
Since At ^ 1, we can apply the Hausdorff-Campbell 
decomposition [6J] via the Zassenhaus formula to lowest 
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order to get 



„iAVoAT iHnAT 



+ 0(At2), (A3) 



where the two orders of the operators are possible since 
they commute with O(At^). The fuU evolution operator 
is then given by 



_ iHi{T-Ti) iAVoAt iHo{To+AT) 



(A4) 



If we take into account that ti = tq + At, wc finally write 

(A5) 



^iHr _ ^iHi{T-Ti) ^IAVqAt ^iHoTi 



If we apply this operator to an initial wave function, we 
see that 



|0(t)) = e''-^i(^-^i)e"^^«'^^e'^«^i|0(O)) 



^iHi(T-Ti) iAVoAr 



\^{ri)). 



(A6) 



It turns out that the presence of an impulse at time ti 
only produces a multiplication by the diagonal operator 
in position space, e''^^"^'''^'^. If we define 



(A7) 



the resulting amplitude can be propagated to future 
times using the harmonic oscillator Feynmann propaga- 
tor in the final medium. 



Appendix B: Symmetry breaking potential 

Close to the origin, |xp = + so we can 

perform a Taylor expansion of the evolution operator in 
(|A7|) in the complex variable x and keep the lower order 
terms. Because of the Cm invariance of the potential, 
there arc only two types of CAr-invariant products of x 
and X* that can appear in this Taylor expansion: xx* — 

Here, we consider discrete rotational symmetry of or- 
der A'^ = 4. If we perform a Taylor expansion on the ar- 
bitrary impulse function V{x) in both variables and keep 
the allowed terms mentioned previously, the potential of 
the impulse can be expanded to read 

AFo(x) - "0 + «i|x|' + "21x1" + v^x^ + v,x*^ + O(x^), 

(Bl) 

where wq, wi, wi, vq and vi are constants. 

This potential presents the most general form of a C4 
invariant potential close to the symmetry axis. Since we 
assume that the first medium is 0(2) invariant, it is clear 



that the only terms that break the symmetry into C4 are 
X^ and x*^- Since we are only analyzing the result of 
the symmetry breaking process, it is sufficient to only 
consider the x^ and x* terms. We take uq = ui = U2 = 
and proceed to evaluate the form of the function after 
the action of the symmetry breaking impulse. 

By only considering the symmetry breaking terms, our 
evolution operator becomes 



giAVoAr _^ ^iArivox'+Vix"^) ^ 



Appendix C: Energy 



(B2) 



Let us evaluate the energy in the 3D formalism. The 
inner product of r' with r" is 

r' -r" ^y-{-^Iy + {y X A)) 

= -nl{v' • r) +r' • (r X A) 
= -Vtl{v' • r) + A - (r' X r) 



= -r22(r' .r) _ A-L. 
We can rearrange for A • L to get 

-A-L = r'-r" + f}2(r'.r). 
The definition of energy is 

E^\{v' ■v') + \nl{vvl 



(CI) 



(C2) 



(C3) 



and then 

dE 1 d 



1 d 



^=2d7(^'-^')+2i(^or-r) 

= i(r".r' + r'.r") + i|:(l^^r.r) 



(C4) 



(C5) 



We use (|C2p and the expression for fig to obtain 



Finally, if we use our value for A • L, we see that the 
change in energy is consistent with our previous analysis 
in (Hal) 



(C7) 



dE 1 A ,2^2 3 A 3 M 1 



dr 



4r3 



4r Vt 



